\defmodule{FaureSequence}

This class implements digital nets or digital sequences formed by the
 first $n = b^k$ points of the Faure sequence in base $b$.
Values of $n$ up to $2^{31}$ are allowed.
% When the base $b$ is not specified, it is taken as the smallest
% prime number greater or equal to the selected dimension.
One has $r = k$.
The generator matrices are
\eq
  \mathbf{C}_j = \mathbf{P}^j\ \mod b
\endeq
for $j=0,\dots,s-1$, where $\mathbf{P}$ is a $k\times k$ upper
 triangular matrix whose entry $(l,c)$ is the number of combinations
of $l$ objects among $c$\latex{, ${c\choose l}$},
for $l\le c$ and is 0 for $l > c$.
% This matrix $\mathbf{P}$ is the transpose of the $k\times k$ \emph{Pascal matrix}.
The matrix $\mathbf{C}_0$ is the identity, $\mathbf{C}_1 = \mathbf{P}$,
and the other $\mathbf{C}_j$'s can be defined recursively via
$\mathbf{C}_j = \mathbf{P} \mathbf{C}_{j-1} \mod b$.
Our implementation uses the recursion
\begin{equation}
\begin{latexonly}
  {c \choose l} = {{c-1} \choose l} + {{c-1} \choose {l-1}}
\end{latexonly}
\begin{htmlonly}
 \mbox{Combination}(c, l)\ = \ \mbox{Combination}(c-1, l) \; + \;
  \mbox{Combination}(c-1, l-1)
\end{htmlonly}
\end{equation}
to evaluate the binomial coefficients in the matrices $\mathbf{C}_j$,
as suggested by Fox\latex{ \cite{rFOX86a} (see also \cite{fGLA04a}, page 301)}.
The entries $x_{j,l,c}$ of $\mathbf{C}_j$ are computed as follows:
\[
\begin{array}{lcll}
 x_{j,c,c} &=& 1             &\quad\mbox{ for } c=0,\dots,k-1,\\[4pt]
 x_{j,0,c} &=& j x_{j,0,c-1} &\quad\mbox{ for } c=1,\dots,k-1, \\[4pt]
 x_{j,l,c} &=& x_{j,l-1,c-1} + j x_{j,l,c-1}
                      &\quad\mbox{ for } 2\le c < l \le k-1, \\[4pt]
 x_{j,l,c} &=& 0      &\quad\mbox{ for } c>l \mbox{ or } l \ge k.
\end{array}
\]

For any integer $m > 0$ and $\nu\ge 0$, if we look at the
vector $(u_{i,j,1},\dots,u_{i,j,m})$ (the first $m$ digits
of coordinate $j$ of the output) when $i$ goes from
$\nu b^m$ to $(\nu+1)b^m - 1$, this vector takes each of its $b^m$
possible values exactly once.
In particular, for $\nu = 0$, $u_{i,j}$ visits each value in the
set $\{0, 1/b^m, 2/b^m, \dots, (b^m-1)/b^m\}$ exactly once, so all
one-dimensional projections of the point set are identical.
However, the values are visited in a different order for
the different values of $j$ (otherwise all coordinates would be identical).
For $j=0$, they are visited in the same
order as in the van der Corput sequence in base $b$.

An important property of Faure nets is that for any integers $m > 0$
and $\nu\ge 0$, the point set
$\{\bu_i$ for $i = \nu b^m,\dots, (\nu+1)b^m -1\}$
is a $(0,m,s)$-net in base $b$.
In particular, for $n = b^k$, the first $n$ points form a
 $(0,k,s)$-net in base $b$.
The Faure nets are also projection-regular and dimension-stationary\latex{
(see \cite{vLEC02a} for definitions of these properties)}.

To obtain digital nets from the \emph{generalized Faure sequence}
\latex{ \cite{rTEZ95a}}, where $\mathbf{P}_j$ is left-multiplied by some
invertible matrix $\mathbf{A}_j$, it suffices to apply an appropriate
matrix scramble (e.g., via \method{leftMatrixScramble}{}).
This changes the order in which $u_{i,j}$ visits its different
values, for each coordinate $j$, but does not change the set of values
that are visited.  The $(0,m,s)$-net property stated above remains valid.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bigskip\hrule\bigskip

\begin{code}
\begin{hide}
/*
 * Class:        FaureSequence
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.hups;\begin{hide}
import umontreal.iro.lecuyer.util.PrintfFormat;

\end{hide}

public class FaureSequence extends DigitalSequence \begin{hide} {

    // Maximum dimension for the case where b is not specified.
    // Can be extended by extending the precomputed array prime[].
    private static final int MAXDIM = 500;

    // For storing the generator matrices for given dim and numPoints.
    private int[][][] v;
\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}

   public FaureSequence (int b, int k, int r, int w, int dim) \begin{hide} {
      init (b, k, r, w, dim);
   }

   private void init (int b, int k, int r, int w, int dim) {
      if (dim < 1)
         throw new IllegalArgumentException
            ("Dimension for FaureSequence must be > 1");
      if ((double)k * Math.log ((double)b) > (31.0 * Math.log (2.0)))
         throw new IllegalArgumentException
            ("Trying to construct a FaureSequence with too many points");
      if (r < k || w < r)
         throw new IllegalArgumentException
            ("One must have k <= r <= w for FaureSequence");
      this.b    = b;
      numCols   = k;
      numRows   = r;
      outDigits = w;
      this.dim  = dim;

      int i, j;
      numPoints = b;
      for (i=1; i<k; i++) numPoints *= b;

      // Compute the normalization factors.
      normFactor = 1.0 / Math.pow ((double) b, (double) outDigits);
//      EpsilonHalf = 0.5*normFactor;
      double invb = 1.0 / b;
      factor = new double[outDigits];
      factor[0] = invb;
      for (j = 1; j < outDigits; j++)
         factor[j] = factor[j-1] * invb;

      genMat = new int[dim * numCols][numRows];
      initGenMat();
   }
\end{hide}
\end{code}
\begin{tabb}
    Constructs a digital net in base $b$,
    with $n = b^k$ points and $w$ output digits,
    in \texttt{dim} dimensions.
    The points are the first $n$ points of the Faure sequence.
    The generator matrices $\mathbf{C}_j$ are $r\times k$.
    Unless, one plans to apply a randomization on more than $k$ digits
    (e.g., a random digital shift for $w > k$ digits, or a linear
    scramble yielding $r > k$ digits), one should
    take $w = r = k$ for better computational efficiency.
    Restrictions: \texttt{dim} $\le 500$ and $b^k \le 2^{31}$.
\end{tabb}
\begin{htmlonly}
   \param{b}{base}
   \param{k}{there will be b^k points}
   \param{r}{number of rows in the generator matrices}
   \param{w}{number of output digits}
   \param{dim}{dimension of the point set}
\end{htmlonly}
\begin{code}

   public FaureSequence (int n, int dim) \begin{hide} {
      if ((dim < 1) || (dim > MAXDIM))
         throw new IllegalArgumentException
            ("Dimension for Faure net must be > 1 and < " + MAXDIM);
      b = getSmallestPrime (dim);
      numCols = (int) Math.ceil (Math.log ((double) n)
                                 / Math.log ((double) b));
      outDigits = (int) Math.floor (Math.log ((double)(1 << (MAXBITS - 1)))
                                 / Math.log ((double)b));
      outDigits = Math.max (outDigits, numCols);
      numRows = outDigits;
      init (b, numCols, numRows, outDigits, dim);
   }
\end{hide}
\end{code}
\begin{tabb}
  Same as \method{FaureSequence}{}\texttt{(b, k, w, w, dim)}
  with base $b$ equal to the smallest prime larger or equal to \texttt{dim},
  and with \emph{at least} \texttt{n} points.
\hrichard{Ce constructeur devrait-il dispara\^\i tre?}
  The values of $k$, $r$, and $w$ are taken as
  $k = \lceil \log_b n\rceil$ and
  $r = w = \max(k, \lfloor 30 / \log_2 b\rfloor)$.
\end{tabb}
\begin{htmlonly}
   \param{n}{minimal number of points}
   \param{dim}{dimension of the point set}
\end{htmlonly}
\begin{code}\begin{hide}

   public String toString() {
      StringBuffer sb = new StringBuffer ("Faure sequence:" +
                  PrintfFormat.NEWLINE);
      sb.append (super.toString());
      return sb.toString();
   }


   public void extendSequence (int k) {
      init (b, k, numRows, outDigits, dim);
   }


   // Fills up the generator matrices in genMat for a Faure sequence.
   // See Glasserman (2004), \cite{fGLA04a}, page 301.
   private void initGenMat() {
      int j, c, l;
      // Initialize C_0 to the identity (for first coordinate).
      for (c = 0; c < numCols; c++) {
         for (l = 0; l < numRows; l++)
            genMat[c][l] = 0;
         genMat[c][c] = 1;
      }
      // Compute C_1, ... ,C_{dim-1}.
      for (j = 1; j < dim; j++) {
         genMat[j*numCols][0] = 1;
         for (c = 1; c < numCols; c++) {
            genMat[j*numCols+c][c] = 1;
            genMat[j*numCols+c][0] = (j * genMat[j*numCols+c-1][0]) % b;
         }
         for (c = 2; c < numCols; c++) {
            for (l = 1; l < c; l++)
               genMat[j*numCols+c][l] = (genMat[j*numCols+c-1][l-1]
                                        + j * genMat[j*numCols+c-1][l]) % b;
         }
         for (c = 0; c < numCols; c++)
            for (l = c+1; l < numRows; l++)
               genMat[j*numCols+c][l] = 0;
      }
/*
      for (j = 0; j < dim; j++) {
     for (l = 0; l < numRows; l++) {
         for (c = 0; c < numCols; c++)
            System.out.print ("  " + genMat[j*numCols+c][l]);
       System.out.println ("");
      }
        System.out.println ("");
  }
*/
   }

/*
   // Fills up the generator matrices in genMat for a Faure net.
   // See Glasserman (2004), \cite{fGLA04a}, page 301.
   protected void initGenMatNet() {
      int j, c, l, start;
      // Initialize C_0 to the reflected identity (for first coordinate).
      for (c = 0; c < numCols; c++) {
         for (l = 0; l < numRows; l++)
            genMat[c][l] = 0;
         genMat[c][numCols-c-1] = 1;
      }
      // Initialize C_1 to the identity (for second coordinate).
      for (c = 0; c < numCols; c++) {
         for (l = 0; l < numRows; l++)
            genMat[numCols+c][l] = 0;
         genMat[numCols+c][c] = 1;
      }
      // Compute C_2, ... ,C_{dim-1}.
      for (j = 2; j < dim; j++) {
         start = j * numCols;
         genMat[start][0] = 1;
         for (c = 1; c < numCols; c++) {
            genMat[start+c][c] = 1;
            genMat[start+c][0] = ((j-1) * genMat[start+c-1][0]) % b;
         }
         for (c = 2; c < numCols; c++) {
            for (l = 1; l < c; l++)
               genMat[start+c][l] = (genMat[start+c-1][l-1]
                                     + (j-1) * genMat[start+c-1][l]) % b;
         }
         for (c = 0; c < numCols; c++)
            for (l = c+1; l < numRows; l++)
               genMat[start+c][l] = 0;
      }
   }
*/

   // Returns the smallest prime larger or equal to d.
   private int getSmallestPrime (int d) {
      return primes[d-1];
   }

   // Gives the appropriate prime base for each dimension.
   // Perhaps should be internal to getPrime, and non-static, to avoid
   // wasting time and memory when this array is not needed ???
   static final int primes[] =
      {2,2,3,5,5,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,
     23,23,23,29,29,29,29,29,29,31,31,37,37,37,37,37,37,41,41,41,
     41,43,43,47,47,47,47,53,53,53,53,53,53,59,59,59,59,59,59,61,
     61,67,67,67,67,67,67,71,71,71,71,73,73,79,79,79,79,79,79,83,
     83,83,83,89,89,89,89,89,89,97,97,97,97,97,97,97,97,101,101,101,
     101,103,103,107,107,107,107,109,109,113,113,113,113,127,127,127,127,127,127,127,
     127,127,127,127,127,127,127,131,131,131,131,137,137,137,137,137,137,139,139,149,
     149,149,149,149,149,149,149,149,149,151,151,157,157,157,157,157,157,163,163,163,
     163,163,163,167,167,167,167,173,173,173,173,173,173,179,179,179,179,179,179,181,
     181,191,191,191,191,191,191,191,191,191,191,193,193,197,197,197,197,199,199,211,
     211,211,211,211,211,211,211,211,211,211,211,223,223,223,223,223,223,223,223,223,
     223,223,223,227,227,227,227,229,229,233,233,233,233,239,239,239,239,239,239,241,
     241,251,251,251,251,251,251,251,251,251,251,257,257,257,257,257,257,263,263,263,
     263,263,263,269,269,269,269,269,269,271,271,277,277,277,277,277,277,281,281,281,
     281,283,283,293,293,293,293,293,293,293,293,293,293,307,307,307,307,307,307,307,
     307,307,307,307,307,307,307,311,311,311,311,313,313,317,317,317,317,331,331,331,
     331,331,331,331,331,331,331,331,331,331,331,337,337,337,337,337,337,347,347,347,
     347,347,347,347,347,347,347,349,349,353,353,353,353,359,359,359,359,359,359,367,
     367,367,367,367,367,367,367,373,373,373,373,373,373,379,379,379,379,379,379,383,
     383,383,383,389,389,389,389,389,389,397,397,397,397,397,397,397,397,401,401,401,
     401,409,409,409,409,409,409,409,409,419,419,419,419,419,419,419,419,419,419,421,
     421,431,431,431,431,431,431,431,431,431,431,433,433,439,439,439,439,439,439,443,
     443,443,443,449,449,449,449,449,449,457,457,457,457,457,457,457,457,461,461,461,
     461,463,463,467,467,467,467,479,479,479,479,479,479,479,479,479,479,479,479,487,
     487,487,487,487,487,487,487,491,491,491,491,499,499,499,499,499,499,499,499,503};

}
\end{hide}
\end{code}
